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Abstract 

In  this  effort  we  investigate  the  behavior  of  a  model  derived  from  homogenization 
theory  as  the  model  solution  in  parameter  estimation  procedures  for  simulated  data  for 
heat  flow  in  a  porous  medium.  We  consider  data  simulated  from  a  model  on  a  perfo¬ 
rated  domain  with  isotropic  flow  and  data  simulated  from  a  model  on  a  homogeneous 
domain  with  anisotropic  flow.  We  consider  both  ordinary  and  generalized  least  squares 
parameter  estimation  procedures. 
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1  Introduction 


Nondestructive  evaluation  is  often  used  to  identify  damage  in  structures,  including  compo¬ 
nents  of  aircraft,  spacecraft,  automobiles,  trains  and  piping,  as  they  age  beyond  their  design 
life  (though  there  are  many  other  uses,  see  [15]).  Many  nondestructive  evaluation  techniques 
(acoustic,  eddy-current,  etc)  detect  and  characterize  damages  through  differences  in  observed 
physical  parameters,  requiring  sophisticated  parameter  estimation  procedures.  In  thermal 
nondestructive  evaluation,  damage  is  often  detected  and  characterized  through  differences 
in  the  observed  thermal  diffusivity.  Nondestructive  evaluation  has  been  developed  in  most 
cases  for  homogeneous  materials.  However  many  structures  of  current  interest  contain  com¬ 
ponents  of  composite  materials  which  are  manufactured  with  a  nontrivial  amount  of  porosity. 
In  [1],  we  considered  modeling  the  flow  of  heat  through  porous  samples  by  use  of  solutions 
of  a  heat  equation  on  a  randomly  perforated  domain.  As  noted,  this  model  was  too  com¬ 
putationally  intensive  for  direct  use  in  parameter  estimation  or  inverse  problem  techniques. 
We  thus  considered  an  approximation  for  a  heat  equation  on  a  perforated  domain  which 
was  derived  through  limiting  processes  in  homogenization  theories  for  the  forward  problem 
in  [2],  Though  the  results  of  [2]  were  encouraging,  good  behavior  of  an  approximation  in 
the  forward  problem  is  not  necessarily  indicative  of  the  behavior  of  the  approximation  in 
inverse  problems.  Moreover,  in  our  earlier  work  we  did  not  consider  random  error  which  is 
often  associated  with  experimental  data.  In  this  paper  we  present  initial  results  suggesting 
that  indeed  the  homogenization  approximate  mathematical  models  developed  in  [2]  and  [1] 
will  perform  well  in  an  inverse  problem  setting.  We  do  this  with  noisy  simulated  data  in 
the  context  of  mathematical  and  statistical  parameter  estimation  procedures  discussed  and 
developed  in  [6],  [13]  and  also  in  [4],  We  treat  data  simulated  with  absolute  error  in  Section  3 
and  data  simulated  with  relative  error  in  Section  4. 


2  Mathematical  Model 

Before  formulating  a  class  of  inverse  problems,  we  consider  several  models  for  the  forward 
problem.  We  first  summarize  a  method  developed  in  [1]  for  modeling  the  flash  heat  exper¬ 
iment  on  a  porous  domain.  We  consider  a  randomly  perforated  domain  Hell,  where  the 
homogeneous,  non  perforated  domain  is  an  L\  x  L2  rectangle  ( Li  is  the  length  in  the  hori¬ 
zontal  direction).  As  convention,  we  take  spatial  coordinates  ( x ,  y),  where  x  is  the  horizontal 
coordinate  and  y  is  the  vertical  coordinate.  We  assume  nr  randomly  placed  pores  Cli  with 
boundaries  dCli  for  i  =  1,  2  . . .  nr,  which  are  generated  using  methods  described  in  [1] ,  [2]  and 
[18].  We  assume  that  these  pores  do  not  intersect  with  each  other  nor  the  boundaries  of  Cl. 
The  perforated  domain  H  is  given  by  H  \  (U’!l1flj).  The  four  boundaries  of  H  which  are  also 
the  four  exterior  boundaries  of  H  are  denoted  u y  for  i  =  1,  2,  3, 4  (as  depicted  in  Figure  2.1). 
We  model  the  flash  heat  experiment  which  approximates  an  experiment  where  the  bottom 
boundary  uq  =  {(x,y) \x  G  [0,Li],y  =  0}  is  heated  by  a  flash  heat  source  [16].  Throughout 
this  document,  we  will  refer  to  uq  as  the  source  boundary.  We  model  the  dynamics  of  the 
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flash  heat  experiment  on  fl  with  the  partial  differential  equation 


'u™nd  -  «V  •  (Wand)  +  7'urand  =  0  in  O  x  (0,  T) 
a^=0  on  dfli  x  (0,  T) 

<adv^  =  _Xu^nd  on  Ud=1  uji  x  (0, T)  (2.1) 

=  S>Z[0,t,]  (*)  -  Awrand  on  a;4  x  (0,  T) 

urand(0,x)  =  0, 


where  a  is  the  thermal  diffusivity  of  the  material  12,  7  corresponds  to  loss  in  the  direction 
orthogonal  to  the  domain  (the  z  direction)  and  A  corresponds  to  loss  on  the  boundary  of 
the  the  rectangle  12. 


0)2 


Figure  2.1:  An  example  randomly  perforated  domain  12  (enlarged  view). 


The  flash  heat  input  is  modeled  by  the  term  SfT\o  ta]  (t)  where 


^[0  ,ts]  (t) 


1,  for  t  <  ts 
0  for  t  >  ts. 


There  are  a  number  of  difficulties  associated  with  using  (2.1)  as  a  model  when  carrying  out 
inverse  problems.  The  computational  time  associated  with  solving  the  forward  problem  (2.1) 
for  -urand  on  a  time  interval  of  interest  in  the  flash  heat  experiments  is  roughly  two  minutes. 
This  is  prohibitively  long  for  use  in  the  inverse  problems  as  well  as  in  some  simulation 
applications  of  the  forward  model.  Beyond  the  computational  intensity  associated  with 
solving  (2.1),  the  random  geometry  of  thin  porous  samples  (which  we  model  as  f2)  is  not 
precisely  known  for  the  nondestructive  evaluation  (NDE)  applications  of  interest  and  thus 
we  cannot  assume  that  is  known  a  priori. 

I11  [2],  we  discussed  the  approximation  of  the  heat  equation  on  the  random  domain 
with  results  derived  from  homogenization  theory  (see  [2], [8],  [11],  [9],  [12]  and  [10]  and  the 
references  therein  for  details).  Using  the  results  of  homogenization  theory,  we  obtain  the 
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limit  system 


'pvUt  -  aV  •  (A°W)  +  ^pvU  =  0  infix  (0,  T) 

<  a^  =  -\U  on  Uf=1  Ui  x  (0,  T) 

'  =  Sfl[0M (t)  -  XU  on  cn4  x  (0,  T) 

U(0,x)  =  0, 


where  py  is  the  proportion  of  0  occupied  by  Cl  (pv  =  area  °(  -s )  and  A0  is  the  2x2  homogenized 

matrix  which  can  be  readily  calculated  using  methods  described  in  [2],  On  the  boundaries 

we  use  the  notation  =  n  ■  A°\7U  where  n  is  the  exterior  unit  normal  vector.  The 

ovA  o 


action  of  A0  is  to  approximate  the  isotropic  flow  through  the  random  domain  fl  around 
the  perforations  Cli  with  anisotropic  flow  through  the  homogeneous  rectangle  hi.  For  more 
discussion  of  this  aspect  of  our  homogenization  approximation,  see  [2],  Specifically,  under 
reasonable  assumptions,  one  has  in  particular  the  following  convergences  as  nr  — >  oo: 


urand  pvU  weakly  in  L2(0;  T;  H\£l)), 
v{urand)  v{U)  weakly  in  L2(0,  T;  L2(o;4)), 


(2.3) 


where  ~  denotes  the  zero  extension  of  a  function  defined  on  to  all  of  Cl,  and  v  is  the  linear 
trace  operator  v  :  L2(0,  T;  if1(Q))  i->  L2(0,  T;  L2(o;4)).  Moreover  (for  more  details  see  [5]),  a 
corrector  result  also  shows  that 

(  lim  || urand  —  f/||c([o,r];L2(n))  =  0, 

I  nr— >oo  (2.4) 

|  lim  \\'Vurand  —  Ca,WU\\c([o,T\-L2(n))  —  0, 

V,  nr^ oo 

where  Cq  is  the  corrector  matrix  associated  to  the  elliptic  problem  corresponding  to  (2.1). 
Actually  one  also  has  error  estimates  of  the  type 

\\urand  -  U\\L2{n)  <  eV2C,  a.  e.  for  t  G  (0,  T). 

Based  on  this  and  on  the  convergences  (2.3),  (2.4),  we  propose  to  use  U,  the  solution  of 
(2.2),  as  a  model  solution  in  the  ordinary  least  squares  estimation  (OLS)  procedures  with 
simulated  data  generated  using  U  with  added  absolute  random  error  and  with  simulated 
data  generated  using  ■urand  (the  solution  of  (2.1))  again  with  added  absolute  random  error  in 
Section  3.  In  Section  4,  we  use  U  as  a  model  solution  in  the  generalized  least  squares  (GLS) 
estimation  procedures  with  simulated  data  generated  using  U  with  added  relative  random 
error  and  with  simulated  data  generated  using  -urand  with  added  relative  random  error. 


3  Ordinary  Least  Squares 

We  assume  a  statistical  model  which  describes  data  with  random  error  that  has  zero  mean,  is 
independent  and  has  constant  variance.  We  discuss  the  OLS  parameter  estimation  procedure 
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in  general  and  then  go  on  to  present  results  from  carrying  out  the  OLS  parameter  estimation 
procedure  on  the  two  different  forms  of  simulated  data  in  Section  3.1. 

We  consider  the  full  parameter  set  9  =  (7,  a,  A),  and  subsets  9X  =  (yA,  a\)  (corresponding 
to  the  assumption  that  the  boundary  loss  parameter  A  is  known)  and  6 7  =  (a7,  A7)  (corre¬ 
sponding  to  the  assumption  that  the  parameter  that  models  loss  in  the  orthogonal  direction 
7  is  known).  Having  introduced  the  three  parameter  sets  of  interest,  we  will  use  without  loss 
of  generality  the  notation  9&  to  represent  any  one  of  the  three  parameter  sets  9,  6X  and  6 7  . 

Full  state  observation  is  rare,  especially  when  the  set  of  states  is  continuous  (for  our 
problem  fsH).  Often,  when  one  performs  thermal  nondestructive  evaluation,  data  is  given 
by  the  output  of  an  IR  camera  on  the  boundary  07.  To  model  the  resulting  pixels,  we  define 
observation  operators 

C<(0)  =  j  f  (f>(s,0)ds. 

J  Xi 

Thus  Ci  yields  the  average  value  of  functions  along  intervals  of  length  l  starting  at  x  =  Xi  on 
07  (the  source  boundary  as  well  as  the  observation  boundary).  We  suppose  the  “perfectly 
resolved”  data  is  given  by 

U,,{Ot)  =  \  r^U(tjts,0-X)ds  (3,1) 

J  Xi 

or  Uijief)  =  Cjl'itj.  •:(){()  where  9 f  is  the  “true”  parameter  value.  We  will  denote  m 
spatial  nodes  Xi  =  Xi,x2,  ■  ■  ■  ,xm,  and  n  temporal  nodes  t3  =  ti,t2,  ■ . .  ,tn.  The  statistical 
assumptions  that  underlie  the  OLS  parameter  estimation  procedure  corresponding  to  this 
observation  process  is  then  that  data  are  given  by  realizations  of  the  random  process  Y  ,3 
which  is  defined  as 

Yij  =  Uij(9*)  +  £ij,  (3.2) 

where  £tJ  is  a  random  variable  that  satisfies  (3.3)  below.  The  random  variable  (random 
error)  £tJ  is  further  assumed  to  have  zero  mean,  be  independent  and  have  constant  variance. 
More  precisely,  we  assume 


E{Sij)  =  0 

Var(£ij )  =  Uq  (3.3) 

Cov(£i:j,£kh )  =  0  for  (i,j)^(k,h). 

It  is  important  to  emphasize  that  Y is  a  random  variable  with  realizations  yt] .  The  real¬ 
ization  yij  would  correspond  to  observed  data. 

In  order  to  obtain  the  parameter  estimate  9#,  we  must  minimize  the  OLS  cost  functional 

m  n 

J(0*)  =  EE  (C«(9#)^B«)2.  (3.4) 

4=1  j=  1 

So  for  each  data  set  {ytj},  the  parameter  estimate  is  given  by 

9 *  =  arg  min  J(9*)  (3.5) 

e#ee# 
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where  0^  is  an  admissible  parameter  set.  The  variance  estimate  is  then  given  by 


■HP*) 

nm  —  p 


(3.6) 


where  p  is  the  number  of  parameters;  thus  p  —  2  for  9 #  =  9 7  and  9 #  =  9X  while  p  —  3  for 
9#  =  9. 

Data  collection  nodes  can  be  selected  in  many  sophisticated  ways  such  as  by  using  SE- 
optimal  design,  E-optimal  design,  D-optimal  and  c-optimal  design  methods  (see  [7]  and  refer¬ 
ences  therein).  However,  this  is  not  the  focus  of  our  current  efforts  so  we  will  simply  examine 

_ip 

the  traditional  sensitivity  functions  to  select  data  collections  nodes.  For  each  parameter  9jk 
in  the  parameter  set  9 #,  the  associated  traditional  sensitivity  function  (the  sensitivity  of  the 
model  solution  with  respect  to  9f)  V9*  =  corresponds  to  the  sensitivity  of  the  solution 

with  respect  to  the  kill  parameter.  In  places  where  the  sensitivity  V9*  is  zero,  one  cannot 
obtain  any  information  about  the  A; t li  parameter.  However,  sensitivity  information  must  be 
used  with  care  in  design  of  inverse  problems  [3].  For  example,  the  information  is  local  in 
nature  (i.e.,  depends  on  the  values  9 #  at  which  the  derivatives  are  evaluated)  and  moreover, 
one  should  not  exclusively  choose  nodes  in  the  regions  of  the  highest  sensitivity. 

We  examined  the  output  sensitivity  functions 


F7(f,oy)  =  (jfx  +  U(t,s,0-(10~3,  2.9167,  0.01))  ds 

Va(t,Xi)  =  ^  Q/  *+  U(t,s,0-,(W~3, 2.9167, 0.01))  ds^j  (3.7) 

FA(f,oy)  =  ^  Q  f  I+  f/(t,s,0;(10-3,  2.9167,  0.01))  ds 

\  J  Xi 


for  t  =  0.57  and  xt  =  0,0.57,1.14,1.71  using  calculations  detailed  in  Appendix  A.  We 
chose  the  values  of  ay  and  t  based  on  an  example  pixel  width.  For  9 #  we  chose  7  =  10~3, 
a  =  2.9167  and  A  =  0.01  based  on  physically  reasonable  values  [2],  These  are  also  the  values 
of  9 #  that  we  used  for  the  simulated  data  in  Sections  3.1  and  4.1. 

After  inspecting  H7,  Va  and  Vx  for  different  values  of  ay,  we  observed  little  difference 
and  thus  we  arbitrarily  chose  two  sets  of  spatial  nodes  ay  G  {0,  0.57}  and  ay  G  (0,  0.57, 1.14} 
to  consider  the  effect  of  sparsity  of  spatial  nodes  on  the  inverse  problem.  As  depicted  in 
Figures  3.1(a)  and  (b),  the  values  of  F7(f,ay),  Va(t,Xi)  and  Vx(t,a y)  vary  over  time.  The 
sensitivity  with  respect  to  a,  Va(t,Xi)  goes  to  zero  quickly  after  the  end  of  flash-heating  at 
ts  =  0.6  s  in  Figure  3.1(a).  There  is  often  more  measurement  error  during  flash-heating  so 
we  take  many  times  just  after  the  end  of  flash  heating  (0.6  +  0.6  +  j|q,  . . .  0.6  +  y|g} 

to  gain  information  about  the  parameter  a.  The  sensitivities  to  7  and  A  are  depicted  in 
Figure  3.1(b)  and  go  to  zero  much  more  slowly  than  Va.  We  take  times  (20,40, . . . ,  140} 
to  gain  information  about  the  parameters  7  and  A.  We  note  that  in  Figure  3.1(b)  the 
sensitivity  with  respect  to  7  is  less  than  the  sensitivity  with  respect  to  A.  This  suggests 
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Sensitivity  of  the  model  solution  to  cv 


(a) 

Sensitivity  of  the  model  solution  to  7  and  A 


(b) 


Figure  3.1:  The  sensitivity  functions  of  jj(t ,  s,  0)  ds  where  U  is  the  solution  of  (2.2) 

with  7  =  10~3,  a  =  2.9167  and  A  =  0.01.  (a)  The  sensitivity  with  respect  to  a  for  time 
t  G  [0,  2],  (b)  The  sensitivity  with  respect  to  A  and  7  for  time  t  G  [0,  300]. 
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that  there  could  be  problems  estimating  7;  we  will  further  discuss  the  well-posedness  of 
estimating  7  in  Section  3.1. 

Using  the  sensitivity  matrix  with  entries 

Xi+m(j-i)Ae)  =  4rm)  (3-8) 

Ot>k  c,=f) 

we  may  use  asymptotic  distribution  results  (see  [6])  to  calculate  the  estimated  covariance 
m.atnX 

no*)  =  n  [xT{0*)x{0*))  (3.9) 

and  the  estimated  standard  error  for  the  kth  parameter 

SE  (§#)  =  yjxkk(6*).  (3.10) 

3.1  Simulated  Ordinary  Least  Squares  Data 

We  would  like  to  consider  data  motivated  by  (2.1)  because  we  suspect  that  this  model  will 
generate  solutions  that  resemble  experimental  data.  We  cannot,  however,  use  (2.1)  as  a 
model  solution  in  the  OLS  cost  functional  (3.4)  because  the  random  geometry  fl  in  (2.1)  is 
not  a  priori  knowledge  in  most  thermal  nondestructive  evaluation  applications.  On  the  other 
hand  we  observe  that  Uij{6 )  of  (2.2)  depends  on  py,  the  porosity  ration  for  which  reasonable 
a  priori  estimates  are  available.  Thus  we  are  concerned  with  the  behavior  of  Uij{6)  as  a 
model  solution  in  the  inverse  problem  with  data  motivated  by  urand.  It  is  not  simple  to 
ensure  that  the  OLS  assumptions  (i.e. ,  model  independent  constant  variance  errors  that  are 
iid  as  in  (3.3))  are  satisfied  even  in  the  most  straight  forward  cases  because  relations  between 
parameters  can  cause  violations  in  the  OLS  assumptions.  For  instance,  the  estimation  of  7  is 
ill-posed.  In  the  partial  differential  equation  (2.2),  if  we  make  the  transformation  U  =  e_7*Z 
in  (2.2)  the  partial  differential  equation  system  becomes 

{pvZt  -  aV  •  (A°VZ)  =0  on  Cl  x  (0,  T) 

011  u?=i  Ui  x  (0,  T )  (3.11) 

aW^  =  ~  AZ  011  x  (0, T), 

where  X[0its](f)  =  e7tZ[0it3](f).  This  means  that  if  we  rewrite  the  cost  functional  in  (3.4)  as 

m  n 

J m  =  J2  E  (cvfo,  •;<*)-  dU  <3-12) 

i=  1  j= 1 

we  can  see  that  the  cost  functional  J{6)  in  (3.12)  where  U(tj,  •;  9)  is  the  solution  of  (2.2)  can 
equivalently  be  written  as  the  weighted  least  squares 

m  n 

•w  =  ££e-27‘'  (CiZ{tj7  •;  9)  -  e^d^)2, 

*=  1  3= 1 
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(3.13) 


and  an  equivalent  OLS  formulation  for  Z  would  drop  the  weight  terms  e^2ltj .  Here  Z(9)  is 
the  solution  of  (3.11)  and  d^  is  the  observed  data.  This  implies  that  by  including  7  in  the 
inverse  problem  formulation,  any  noise  in  the  data  might  be  amplified.  This  will  also  cause 
model  dependence  in  the  error  in  d,7  which  violates  the  error  assumptions  of  (3.3)  in  the 
corresponding  OLS  formulation.  This  is  evidence  that  the  parameter  7  plays  a  very  different 
role  in  the  estimation  problems  when  compared  to  that  of  either  a  or  A,  and  along  with  the 
sensitivity  results  of  Figure  3.1,  suggests  an  additional  level  of  ill-posedness  in  estimating  7 
that  is  not  present  in  problems  where  7  is  held  fixed. 

Because  of  the  subtlety  in  verifying  the  OLS  error  assumptions  (3.3)  directly,  we  will 
compare  data  simulated  using  the  solution  of  (2.2)  with  data  simulated  using  the  solution  of 
(2.1).  We  simulate  data  motivated  by  OLS  assumptions  and  the  solution  of  (2.2)  with  the 
random  process 

D i:j(a)  =  Uij{  1(T3,  2.9167,  0.01)  +  oB^  (3.14) 

where  Br]  is  a  random  variable  which  follows  a  standard  normal  distribution  or  B{j  ~ 
0,  l2).  We  consider  two  sets  of  spatial  nodes  27  =  0,  0.57  and  27  =  0,  0.57, 1.14. 

We  simulate  data  motivated  by  the  OLS  error  assumptions  (3.3)  and  the  solution  of  (2.1) 
with  the  random  process 

D  ™nd(a)=n^d  +  aB,ij  (3.15) 

where  u™nd  is  given  by 


rand  _  _ 

L  0 


rxi+i 


uraXLd(tj ,  s,  0;  (10~3,  2.9167,  0.01))  ds, 


'Xi 


(3.16) 


with  uiand(tj,s,  0;  (10~3,  2.9167,  0.01))  the  solution  of  (2.1)  with  (7,  a,  A)  =  (10“3,  2.9167,  0.01) 
on  a  randomly  perforated  geometry  O. 

For  realizations  of  D ij(cr)  and  D^nd(cr)  each  data  set  is  analyzed  using  the  results  of  OLS 
asymptotic  theory  for  the  parameter  sets  9 7  =  (a7,  A7),  9X  =  (7^,  07)  and  6  =  (7,  a,  A).  So  for 
each  data  set  we  preform  three  inverse  problems  calculating  three  parameter  estimates  d7,  9X 
and  9  using  (3.5),  three  variance  estimates  b2,  a\  and  cr2  using  (3.6),  three  covariance  matrix 
estimates  E(07),  S(dA),  and  E(0)  using  (3.9),  and  the  standard  error  for  each  parameter  in 
the  sets  #7,  9X  and  9,  denoted  SE(d7),  SE(A7),  SE(7A),  SE(dA),  and  SE(7),  SE(d),  SE(A) 
using  (3.10). 

There  are  many  different  ways  to  consider  these  parameter  estimation  and  uncertainty 
quantification  problems.  We  will  consider  the  difference  between  the  parameter  estimate  and 
the  “true”  parameter  values  70  =  10”3,  ao  =  2.9167,  and  Ao  =  0.01  which  will  give  us  insight 
into  the  accuracy  of  the  parameter  estimate  Of.  Recall  the  ratio  SE 0f)/Of  is  related  to 

^-LL 

the  uncertainty  associated  with  the  parameter  estimate  9J^ .  Thus  we  also  consider  the  ratio 
of  the  estimated  standard  error  to  the  parameter  estimate.  When  this  ratio  is  large,  there  is 
little  confidence  in  the  value  of  the  parameter  §f.  For  instance,  the  95%  confidence  interval 
[6]  for  rn  —  3  for  7  is  given  by  (7  —  2.02  x  SE(7),7  +  2.02  x  SE(7)).  In  this  example  when 
SE(7)/7  is  greater  than  0.5,  the  confidence  interval  actually  covers  possible  negative  values 
for  7. 
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We  compare  the  results  of  ordinary  least  squares  parameter  estimation  procedure  for 
realizations  of  Y),}](a)  with  the  results  of  ordinary  least  squares  parameter  estimation  proce¬ 
dure  for  realizations  of  D™nd(u).  In  Figures  3. 2-3. 7,  the  simulated  data  is  taken  at  temporal 
collection  nodes  tj  =  0.6  +  y^,  0.6  +  y|g, . . . ,  0.6  +  20, 40, ... ,  140  and  spatial  nodes 

Xi  =  0,  0.57, 1.14.  Figure  3.2(a)  depicts  the  difference  a#  —  ao  f°r  each  of  five  realizations  at 
each  level  of  added  noise  a  =  0.015,  0.030, . . . ,  0.090  for  (<r).  The  estimates  for  the  three 
parameter  subsets  do  not  appear  to  vary  much.  That  is,  we  see  that  for  each  realization  the 
differences  a  —  a0,  a7— a0  and  a\  — a0  are  relatively  close.  Figure  3.2(b)  depicts  the  difference 
d#  —  cto  for  each  of  five  realizations  of  D))md  (a)  at  each  level  of  noise  o  =  0,  0.015, . . . ,  0.090, 
we  see  that  the  differences  d#  —  ao  again  do  not  appear  to  vary  much  for  each  realization. 
Moreover,  the  results  depicted  in  Figure  3.2(a)  resemble  the  results  depicted  in  Figure  3.2(b). 
This  suggests  that  using  realizations  of  D^ct)  versus  D)“ld((r)  does  not  affect  the  accuracy 
of  the  estimates  a#.  In  Figures  3.3  (a)  and  (b),  we  examine  the  ratios  SE(d#)/d#  for  the 
five  realizations  of  T>l3(a)  and  D^nd(cr),  respectively.  In  Figure  3.3(a),  we  see  that  the  un¬ 
certainty  associated  with  the  estimate  d  is  larger  than  the  uncertainty  associated  with  the 
estimates  d7  and  a\.  It  also  appears  that  the  ratio  SE(d#)/d#  varies  linearly  with  a.  These 
observations  are  valid  for  the  ratios  depicted  in  Figure  3.3(b),  as  well.  The  similarities  in 
Figure  3.3(a)  and  Figure  3.3(b)  suggest  that  there  is  not  a  significant  difference  in  using 
realizations  Dq(o-)  and  DUnd(cr)  regarding  the  uncertainty  associated  the  OLS  parameter 
estimate  a#. 

We  also  considered  the  parameter  A.  Figures  3.4(a)  are  the  differences  A#— A0  (A0  —  0.01) 
for  the  realizations  used  in  Figures  3.2(a)  and  3.3(a).  Unlike  the  parameter  a,  there  is  a 
large  difference  in  the  A#  —  Ao-  The  magnitude  of  the  difference  A  —  Ao  is  much  larger  than 
the  magnitude  of  the  difference  A7  —  A0  indicating  that  estimating  7  adds  inaccuracy  to  the 
estimate  of  A.  In  Figure  3.4(b),  we  have  plotted  the  difference  A#  —  Ao  for  the  realizations 
depicted  in  Figures  3.2(b)  and  3.3(b).  In  Figure  3.4(b),  we  see  that  the  differences  A7  —  Ao 
and  A  —  Ao  for  the  realizations  of  D^nd(cr)  resemble  these  differences  for  the  realizations  of 
D,  j(a)  in  Figure  3.4(a)  so  we  suspect  that  the  error  associated  with  the  approximation  of 
■urand  j'ypg  soiution  of  (2.1))  with  U  (the  solution  of  (2.2))  in  the  model  solution  does  not 
affect  the  estimate  of  A7  nor  that  of  A. 

In  Figure  3.5(a),  we  see  that  the  ratio  SE(A7)/A7  appears  to  be  linear  in  a  for  the 
realizations  of  D ij(a).  This  linear  pattern  is  similar  to  the  linearity  in  Figure  3.5(c)  which 
depicts  the  ratio  SE(A7)/A7  for  the  realizations  of  D™nd(cr).  Moreover,  the  ratio  SE(A7)/A7 

is  on  about  the  same  scale  in  Figures  3.5(a)  and  (c).  So  for  A7,  there  does  not  appear  to  be  a 
significant  difference  in  using  data  generated  by  (3.14)  versus  data  generated  by  (3.15)  in  the 
accuracy  of  the  parameter  estimate  A7  nor  in  the  uncertainty  associated  with  the  parameter 
estimate  A7. 

In  order  to  consider  the  ratio  SE(A)/A,  we  plotted  the  logarithm  of  this  quantity  in 
Figures  3.5(b)  (for  realizations  of  Djj(cr))  and  3.5(d)  (for  realizations  of  D)and(cr))  because 
in  both  examples  these  quantities  vary  greatly.  The  ratio  SE(A)/A  in  both  Figure  3.5(b) 
and  Figure  3.5(d)  appears  to  grow  exponentially  with  added  error  cr  which  is  not  surprising 
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Figure  3.2:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  y^q)  0.6  +  y§q,  •  •  • ,  Y20  ’  40, . . . ,  140  and  spatial  nodes  x,  =  0,  0.57, 1.14.  The  points 

denoted  with  0  are  the  difference  d7  —  ao,  the  points  denoted  are  the  a\  —  ao,  the  points 
denoted  *  are  the  difference  a  —  ao  (a)  The  result  of  five  realizations  of  D^(cr)  for  values 
of  cr  =  0.015,  0.030, ...,  0.090  (b)  The  result  of  five  realizations  of  D™nd(c)  for  values  of 
<7  =  0,  0.015, . . . ,  0.090. 
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Figure  3.3:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  y^q;  0.6  +  y§q,  •  •  • ,  ^20’  20, 40, . . . ,  140  and  spatial  nodes  x,  =  0,  0.57, 1.14.  The  points 
denoted  with  0  are  the  ratio  SE(d7)/d7,  the  points  denoted  are  the  SE(dA)/dA,  the  points 
denoted  *  are  the  difference  SE(d)/d  (a)  The  result  of  five  realizations  of  D,7(ct)  for  values 
of  cr  =  0.015,  0.030, ...,  0.090  (b)  The  result  of  five  realizations  of  D™nd(c)  for  values  of 
<7  =  0,  0.015, . . . ,  0.090. 
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based  on  our  observation  in  (3.13)  that  7  in  the  parameter  estimation  increases  the  error  in 
the  data. 

The  estimation  of  7  is  ill-posed  for  the  parameter  values,  spatial  nodes  and  temporal  nodes 
in  this  problem.  In  Figures  3.6(a)  and  (b),  we  see  that  the  difference  7  —  70  (70  =  1CT3)  is 
very  large  compared  to  70.  This  is  true  for  both  the  realizations  of  D ij(cr)  in  Figure  3.6(a) 
and  the  realizations  of  D)|nd(a)  in  Figure  3.6(b).  Though  the  differences  7a  —70  are  smaller 
than  7  —  70  Figures  3.6(a)  and  (b),  in  Figures  3.7(a)  and  (c)  we  see  that  the  uncertainty 
associated  with  the  estimate  7a  is  ver  large.  In  Figure  3.7(a),  we  see  that  for  the  realizations 
of  D,  j(a)  the  ratios  SE(7a)/7a  does  not  appear  to  grow  exponentially  with  a  but  the  values 
of  SE(7a)/7a  are  too  large  to  be  confident  of  even  the  sign  of  The  use  of  D-“ld(a)  in 
the  OLS  parameter  estimation  procedure  does  appear  to  affect  the  uncertainty  associated 
with  the  parameter  estimate  7a  as  we  see  in  the  exponential  growth  of  SE(7a)/7a  with  a  for 
realizations  of  D™nd (a)  in  Figure  3.7(c).  The  ratio  SE(7)/7  varies  on  an  exponential  scale 
for  both  realizations  of  Djj(cr)  and  D)|nd(a)  in  Figures  3.7(b)  and  (d),  respectively. 

We  also  considered  data  at  temporal  nodes  tj  =  0.6  +  y^,  0.6  +  . . . ,  0.6  +  ^,20, 

40, . . . ,  140  and  spatial  nodes  xt  =  0,  0.57  to  understand  how  sparsity  of  spatial  data  affects 
the  inverse  problem.  For  each  level  of  added  noise  cr  =  0.015,  0.030, . . . ,  0.090,  we  simulated 
five  realizations  of  Dij(a)  using  (3.14).  We  also  simulated  five  realizations  of  d™nd(<r)  using 
(3.15)  for  each  level  of  added  noise  a  —  0,  0.015, . . . ,  0.090.  For  each  realization  of  D)|nd(<r) 
a  different  random  geometry  is  used  to  solve  (2.1)  for  in  (3.16).  Results  from  these 
simulations  are  given  in  Figures  3.8-3.13. 

For  the  realizations  of  Djj(cr),  the  differences  d#  —  ao  are  plotted  in  Figure  3.8(a)  while 
these  differences  for  the  realizations  of  D™nd(<r)  are  plotted  in  Figure  3.8(b).  Much  like  these 
quantities  in  Figures  3.2(a)  and  (b),  for  each  realization  the  differences  07  —  a o>  d7  —  a0 
and  a  —  07  remain  relatively  close  and  are  on  about  the  same  scale  as  in  Figures  3.2(a) 
and  (b).  The  ratios  SE(d#)/d#  for  the  realizations  of  Dy- (cr)  and  D^nd(cr)  are  plotted  in 
Figures  3.9(a)  and  (b),  respectively.  In  both  figures,  it  appears  that  the  ratio  SE(d)/d 
is  larger  than  the  ratios  SE(d7)/d7  and  SE(b7)/b7.  It  also  appears  that  the  relationship 
between  the  ratios  SE(d#)/d#  and  a  is  linear  as  we  observed  in  Figures  3.3(a)  and  (b)  as  well. 
This  suggests  that  for  the  temporal  and  spatial  nodes  that  we  are  considering  estimating  a 
is  well-posed  and  there  is  little  difference  between  using  realizations  of  D)|nd(a)  and  D ij(cr) 
in  the  estimation  of  the  parameter  a  and  the  estimation  of  the  uncertainty  associated  with 
a. 

In  Figures  3.10(a)  and  (b),  we  see  that  the  difference  A  —  A0  is  larger  than  A7  —  A0  for  the 
realizations  of  D^cr)  and  D)|nd(cr),  respectively.  Figure  3.10(b)  resembles  Figure  3.10(a) 
which  suggests  that  using  realizations  of  D)|nd(a)  rather  than  realizations  of  D ij(cr)  does 
not  have  a  substantial  effect  on  the  accuracy  of  the  parameter  estimates  A#.  The  ratios 
SE(A7)/A7  appear  to  be  similar  for  realizations  of  D p-(cr)  (in  Figure  3.11(a))  and  realizations 
of  D™nd  (a)  (in  Figure  3.11(c))  and  seem  to  vary  linearly  with  a.  The  ratios  SE(A)/A  are 
much  larger  than  the  ratios  SE(A7)/A7,  so  we  plotted  log(SE(A)/A)  for  the  realizations  of 
D ij(a)  in  Figure  3.11(b)  and  for  the  realizations  of  D)|nd((r)  in  Figure  3.11(d).  This  implies 
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Figure  3.4:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT)  ’  0.6  +  ifo  ’ ' ' '  ’  i!o  ’  40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57, 1.14.  The  points  denoted 

with  °  are  the  difference  A7  —  Ao,  and  the  points  denoted  *  are  the  difference  A  —  Ao  (a)  The 
result  of  five  realizations  of  D^- (cr)  for  values  of  cr  =  0.015,  0.030, . . . ,  0.090  (b)  The  result  of 
five  realizations  of  D-|nd(a)  for  values  of  a  —  0,  0.015, . . . ,  0.090. 


14 


SE(A7)/A7  SE(A7)/i 


Ratio  SE(A7)/A7  for  x^  =  0,0.57,1.14 


Log  of  the  ratio  SE(A)/A  for  Xi  =  0,0.57, 1.14 


0.05 

0.045 

0.04 

0.035 

0.03 

0.025 

0.02 

0.015 

0.01 


* 

* 


t 

I 


+ 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 


W 

Ratio  SE(A7)/A7  for  xt  =  0,0.57,1.14 


0.05 


0.04 


0.03 


0.02 


0.01 


* 


+ 


t 

# 


* 


ot_ 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 

a 

(c) 


8 


7- 


'  5  - 


4*  i  +  4 

t  I  m 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 


(b) 


Log  of  the  ratio  SE(A)/A  for  Xi  =  0,0.57, 1.14 


W  2- 


* 

+ 


4 


+ 


4 


+ 


4 

i  i _ i _ i _ i _ i _ i _ i _ i _ i _ j_ 

0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09 


(d) 


Figure  3.5:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  0.6  +  y§q ,  •  •  • ,  Y20 ,  40, . . . ,  140  and  spatial  nodes  =  0,  0.57, 1.14.  (a)  The  ratio 

SE(A7)/A7  for  five  realizations  of  Djj(cr)  for  values  of  cr  =  0.015,  0.030, . . . ,  0.090  (b)  The  log 
of  the  ratio  log  SE(A)/A  for  five  realizations  of  Dy  (cr)  for  values  of  a  —  0.015,  0.030, . . . ,  0.090. 
(c)  The  ratio  SE(A7)/A7  for  five  realizations  of  D^nd(cr)  for  values  of  a  =  0,  0.015, . . . ,  0.090. 


(d)  The  log  of  the  ratio  logSE(A)/A  for  five  realizations  of  D-®'nd(cr)  for  values  of  cr  = 
0,0.015,  0.030,...,  0.090. 
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Figure  3.6:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT) !  0.6  +  ifo  ’  •  •  •  >  120  ’  ^0;  40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57, 1.14.  The  points  denoted 
with  °  are  the  difference  7a  —  7o,  and  the  points  denoted  *  are  the  difference  7  —  70  (a)  The 
result  of  five  realizations  of  D ij(cr)  for  values  of  cr  =  0.015,  0.030, . . . ,  0.090  (b)  The  result  of 
five  realizations  of  D^nd(cr)  for  values  of  a  =  0,  0.015, . . . ,  0.090. 
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Figure  3.7:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  y^Q,  0.6  +  y§q ,  ■  ■  ■ ,  ^20 j  20, 40, . . . ,  140  and  spatial  nodes  x j  =  0,  0.57, 1.14.  (a)  The  ratio 
SE(7a)/7a  for  five  realizations  of  D^(cr)  for  values  of  o  =  0.015,  0.030, . . . ,  0.090  (b)  The  log 
of  the  ratio  log  SE(7)/7  for  five  realizations  of  D?:j  (a)  for  values  of  cr  =  0.015,  0.030, . . . ,  0.090. 
(c)  The  log  of  the  ratio  log(SE(7A)/7A)  for  five  realizations  of  D-^nd(cr)  for  values  of 
a  —  0,  0.015, . . . ,  0.090.  (d)  The  log  of  the  ratio  logSE(7)/7  for  five  realizations  of  D-ynd(a) 
for  values  of  a  =  0,  0.015,  0.030, . . . ,  0.090. 
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a#  —  a o  for  Xi  =  0, 0.57 


Figure  3.8:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  jdg,  0.6  +  ^§0,  ■  ■  ■ ,  Y20’  20, 40, . . . ,  140  and  spatial  nodes  x,  =  0,  0.57.  The  points  de¬ 
noted  with  °  are  the  difference  b7  —  ao,  the  points  denoted  are  the  6>\  —  ao,  the  points 
denoted  *  are  the  difference  a  —  ao  (a)  The  result  of  five  realizations  of  Djj(cr)  for  values 
of  cr  =  0.015,  0.030, ...,  0.090  (b)  The  result  of  five  realizations  of  D  ■|nd(cr)  for  values  of 
<7  =  0,  0.015, . . . ,  0.090. 
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Figure  3.9:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  0.6  +  y§q,  ■  ■  ■ ,  yIo’  20, 40, . . . ,  140  and  spatial  nodes  x,  =  0,  0.57.  The  points  de¬ 

noted  with  0  are  the  ratio  SE(d7)/d7,  the  points  denoted  are  the  SE(dA)/dA,  the  points 
denoted  *  are  the  difference  SE(d)/d  (a)  The  result  of  five  realizations  of  D?7(ct)  for  values 
of  cr  =  0.015,  0.030, ...,  0.090  (b)  The  result  of  five  realizations  of  D^nd(c)  for  values  of 
<7  =  0,  0.015, . . . ,  0.090. 
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that  for  realizations  of  Dy(cr)  and  D-“ld(<j)  estimating  7  causes  a  dramatic  increase  in  the 

uncertainty  associated  with  the  parameter  estimate  A  for  spatial  nodes  ay  =  0,  0.57.  Note 
that  this  effect  appears  to  be  similar  for  realizations  of  Djj(cr)  in  Figure  3.11(b)  as  for 
realizations  of  D™nd (a)  in  Figure  3.11(d). 

The  estimation  of  the  parameter  7  appears  to  be  ill-posed  for  the  spatial  nodes  07  = 
0,0.57,  especially  when  estimating  the  parameter  set  (7,0,  A).  The  differences  7  —  70  are 
several  orders  of  magnitude  larger  than  the  “true”  parameter  70  =  10~3  in  Figure  3.12(a) 
(for  realizations  of  Dy(ff))  and  in  Figure  3.6(b)  (for  realizations  of  D  rv))- 

For  the  spatial  nodes  07  =  0,  0.57,  the  uncertainty  associated  with  the  parameter  es¬ 
timate  7#  is  very  large.  For  all  of  the  examples  of  SE(7#)/7#  in  Figures  3.13(a)-(d)  we 
plotted  log(SE(7#)/7#)  because  the  variation  of  SE(7^)/y#  was  so  large  for  every  exam¬ 
ple.  Figures  3.13(a)  and  (c)  depict  SE(7a)/7a  for  the  realizations  of  (a)  and  D™nd(a), 
respectively.  The  ratios  SE(7a)/7a  are  on  an  exponential  scale  for  realizations  of  Dy(cr)  in 
Figure  3.13(a)  with  spatial  nodes  07  =  0,  0.57  while  the  ratios  SE(7a)/7a  are  on  a  linear  scale 
for  realizations  of  D y(cr)  in  Figure  3.7(a)  which  indicates  that  sparsity  of  spatial  collection 
nodes  affects  the  uncertainty  associated  with  the  parameter  estimates  7#. 
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Figure  3.10:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT)  ’  0-6  +  lio  >  ■  ■  ■  >  i2o>  20, 40, . . . ,  140  and  spatial  nodes  x*  =  0,  0.57.  The  points  denoted  with 
°  are  the  difference  A7  —  Ao,  and  the  points  denoted  *  are  the  difference  A  —  Ao  (a)  The  result 
of  five  realizations  of  D ij(cr)  for  values  of  a  —  0.015,  0.030, . . . ,  0.090  (b)  The  result  of  five 
realizations  of  D™nd(a)  for  values  of  a  =  0,  0.015, . . . ,  0.090. 
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Figure  3.11:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
lio’  0.6  +  y§q 7  •  •  ■  j  120’  20, 40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57.  (a)  The  ratio  SE(A7)/A7 
for  five  realizations  of  Djj(cr)  for  values  of  cr  =  0.015,  0.030, ...,  0.090  (b)  The  log  of  the 
ratio  logSE(A)/A  for  five  realizations  of  Dp- (a)  for  values  of  cr  =  0.015,  0.030, . . . ,  0.090.  (c) 
The  ratio  SE(A7)/A7  for  five  realizations  of  D'J nd (a)  for  values  of  cr  =  0,0.015, . . .  ,0.090. 
(d)  The  log  of  the  ratio  logSE(A)/A  for  five  realizations  of  D.|nd(cr)  for  values  of  cr  = 
0,0.015,  0.030,...,  0.090. 
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Figure  3.12:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
i^o  ’  0-6  +  ifo  >  ■  ■  ■  j  lio  ’  40, . . . ,  140  and  spatial  nodes  27  =  0,  0.57.  The  points  denoted  with 

°  are  the  difference  7a  —  7o,  and  the  points  denoted  *  are  the  difference  7  —  70  (a)  The  result 
of  five  realizations  of  Dp- (a)  for  values  of  a  =  0.015,  0.030, . . . ,  0.090  (b)  The  result  of  five 


realizations  of  D^nd(a)  for  values  of  a  =  0,  0.015, . . . ,  0.090. 
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Figure  3.13:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT)’  b-6  +  y§q>  •  ■  ■ ,  720 ’  20, 40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57.  (a)  The  ratio  SE(7A)/7a 
for  five  realizations  of  D^(cr)  for  values  of  a  =  0.015,  0.030, ...,  0.090  (b)  The  log  of  the 
ratio  logSE(7)/7  for  five  realizations  of  Dy(cr)  for  values  of  a  =  0.015,  0.030, ...,  0.090. 
(c)  The  log  of  the  ratio  log(SE(7A)/7A)  for  five  realizations  of  D™nd(cr)  for  values  of  a  = 
0,  0.015, . . . ,  0.090.  (d)  The  log  of  the  ratio  log  SE(7)/7  for  five  realizations  of  D™nd(<r)  for 
values  of  a  =  0,  0.015,  0.030, . . . ,  0.090. 
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4  Generalized  Least  Squares 


The  generalized  least  squares  (GLS)  parameter  estimation  procedure  like  the  OLS  parameter 
estimation  procedure  is  based  on  an  underlying  statistical  model.  The  error  for  GLS  is 
assumed  to  be  relative  or  proportional  to  the  model  value.  Observations  are  assumed  to  be 
realizations  of  the  random  process  Y i3  given  by 

Ylj  =  Ulj(0  „#)(l  +  £«),  (4.1) 


where  is  assumed  to  have  constant  variance,  zero  mean,  and  mutually  independent  or 


E(Si3) 

Var{£ij ) 
Cov(£ij,  £kh) 


0 

0  for  (i,j)  ±  ( k,h ). 


(4.2) 


Note  that  in  (4.1),  the  error  is  given  by  Uij{df)£i3  so  it  is  proportional  to  the  model  and 
the  variance  is  proportional  to  Assuming  the  statistical  model  given  by  (4.1),  the 

GLS  parameter  estimation  procedure  involves  the  minimization  of  the  cost  functional 


f  yvy  -  \ 2 

The  GLS  parameter  estimate  is  given  by 

9 #  =  arg  min  J(9#), 
e#&e#  v  ' 


(4.3) 


(4.4) 


where  J(6 #)  is  defined  in  (4.3).  We  used  the  iteratively  reweighted  least  squares  method  as 
described  in  [6]  and  [13]  to  minimize  (4.3).  The  GLS  variance  estimate  is  given  by 


j(P*) 

nm  —  p 


(4.5) 


where  again  J{9 #)  is  defined  in  (4.3).  The  nm  x  nm  matrix  of  weights  W(6 #)  has  entries 

1 


1)  (9  ) 


U?A9#) 


(4.6) 


for  i  —  1,  2, . . . ,  m  and  j  —  1,2, ,  n.  The  GLS  covariance  matrix  estimate  S(d^)  is  given 
[6]  by 

E{0*)  =  n  (xT{6*)W(9*)x(0*))  (4.7) 

where  %($#)  is  the  matrix  of  sensitivities  with  entries  given  in  (3.8).  The  standard  error 
estimates  are  again  given  by  the  square  roots  of  the  diagonal  entries  of  the  covariance  matrix 


SE(9f)  =  \JEJF)- 


(4.8) 
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4.1  Simulated  Generalized  Least  Squares  Data 

As  in  the  previous  section,  we  investigate  the  well-posedness  of  the  GLS  parameter  estimation 
procedure  by  considering  the  resulting  parameter  estimates  and  uncertainty  estimates  for 
simulated  data.  We  simulate  data  motivated  by  the  GLS  error  assumptions  in  (4.1)  and 
(4.2)  using  both  the  solutions  ■urand  of  (2.1)  and  solutions  U  of  (2.2).  In  order  to  consider  U 
the  solution  of  (2.2)  as  a  model  solution  in  the  inverse  problem,  we  will  simulate  data  which 
is  given  by  realizations  of  Dp  where  Dp  is  given  by 

D ij(a)  =  Uij {10~3,  2.9167,  0.01)  (1  +  n  B,, )  (4.9) 

where  £/p(10-3,  2.9167,  0.01)  is  given  by  (3.1),  and  I?p  follows  a  standard  normal  distribution 
or  B  ij  rs_/  J\f( 0,  l2).  Again,  as  in  Section  3.1,  we  consider  spatial  nodes  x%  =  0,0.57, 1.14  and 
Xi  =  0,  0.57  and  temporal  nodes  tj  =  0.6  +  y^,  0.6  +  y|p, . . . ,  0.6  +  20, 40, . . . ,  140.  We 

analyzed  five  realizations  of  D p(cr)  for  each  value  of  a  =  0.02,  0.05,  0.10. 

We  also  consider  data  which  is  generated  using  solutions  of  (2.1)  and  motivated  by  the 
GLS  assumptions  with  realizations  of  the  random  process 

D  ™d(a)=u™d(l+aBij),  (4.10) 

where  u™nd  is  defined  in  (3.16),  and  is  a  random  variable  sampled  from  a  standard 
normal  distribution  or  Bij  ~  J\f( 0,  l2).  We  calculated  five  realizations  of  Dpnd(<r)  for  each 
value  of  a  —  0,  0.02,  0.05,  0.10. 

For  both  sets  of  simulations,  we  calculate  the  parameter  estimates  or  6X  =  (77,07), 
01  =  (d7,  A7)  and  6  =  (7,0,  A)  as  defined  in  (4.4)  by  minimizing  (4.3).  We  also  calculate 
SE(d#)  using  (4.8). 

As  in  Section  3.1,  we  consider  the  accuracy  of  the  inverse  problem  by  investigating 
(9#  —  Of  and  the  uncertainty  associated  with  the  inverse  problem  by  investigating  the 
ratios  SE ($#)/$#.  In  Figures  4. 1-4.6,  we  report  on  these  simulations  with  spatial  nodes 
Xi  =  0,  0.57, 1.14.  Figure  4.1(a)  depicts  a#  —  «o  for  five  realizations  of  Dp(<r)  for  each  value 
of  a  =  0.02,  0.05,  0.10.  The  values  of  d#  —  «o  hr  Figure  4.1(b)  for  Eve  realizations  of  D™d(a) 
for  each  value  of  o  =  0,  0.02,  0.05,  0.10  appear  to  be  smaller  than  those  in  Figure  4.1(a).  This 
suggests  that  the  GLS  parameter  estimation  procedure  predicts  a#  more  accurately  for  data 
generated  by  realizations  D)|nd(o-)  than  for  data  generated  by  realizations  of  Dp(cr). 

To  consider  the  uncertainty  associated  with  the  parameter  estimates  a #  —  a0,  we  plotted 
the  ratio  SE(d#)/d#  in  Figure  4.2(a)  for  realizations  of  Dp(cr)  and  in  Figure  4.2(b)  for 
realizations  of  Dpnd(cr).  The  ratios  SE (d#)/d#  appear  to  be  linear  in  a  for  realizations  of 
Dp  (a)  and  D™d(ff)  with  similar  slopes.  It  appears  that  there  is  little  difference  between  us¬ 
ing  realizations  of  Dp(cr)  and  realizations  of  Dpnd(<r)  in  the  GLS  estimate  of  the  uncertainty 
associated  with  the  parameter  estimate  a#. 

We  plotted  the  differences  A#  —  Ao  for  the  realizations  of  Dp(cr)  in  Figure  4.3(a)  and  the 
realizations  of  Dpnd(cr)  in  Figure  4.3(b).  The  differences  between  A7  and  Ao  are  very  small  for 
both  realizations  of  Dp(cr)  in  Figure  4.3(a)  and  realizations  of  D™d((j)  in  Figure  4.3(b).  The 
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a#  —  ao  for  x,  =  0, 0.57, 1.14 
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Figure  4.1:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT) !  0.6  +  ifo  >  •  •  •  >  120  ’  ^0;  40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57, 1.14.  The  points  denoted 
with  °  are  the  difference  d7  —  ao,  the  points  denoted  are  the  a\  —  ao,  the  points  denoted 
*  are  the  difference  d  —  ao  (a)  The  result  of  five  realizations  of  for  values  of  a  — 

0.02,  0.05,  0.10  (b)  The  result  of  five  realizations  of  D^nd(cr)  for  values  of  a  =  0,  0.02,  0.10. 
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SE(d#)/d#  for  Xi  =  0,  0.57, 1.14 
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Figure  4.2:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  y^q)  0.6  +  y§q,  •  •  • ,  Y20  ’  ■  ■  ■ ,  140  and  spatial  nodes  x,  =  0,  0.57, 1.14.  The  points 

denoted  with  0  are  the  ratio  SE(d7)/d7,  the  points  denoted  are  the  SE(dA)/dA,  the  points 
denoted  *  are  the  difference  SE(d)/d  (a)  The  result  of  five  realizations  of  D^(cr)  for  val¬ 
ues  of  cr  =  0.02,0.05,0.10  (b)  The  result  of  five  realizations  of  D™nd(cr)  for  values  of 
a  =  0,0.02,0.05,0.10. 


differences  A  —  A0  (depicted  in  Figure  4.3(a)  for  realizations  of  Dp(cr)  and  Figure  4.3(b)  for 
realizations  of  D)|nd(<j))  have  much  larger  magnitudes  than  the  magnitudes  of  the  differences 
A7  —  Ao-  This  suggests,  as  discussed  in  Section  3.1  (see  the  discussions  involving  (3.11)- 
(3.13)),  that  estimating  7  detracts  from  the  accuracy  of  the  GLS  estimate  of  the  parameter 
A. 

We  further  see  the  effect  of  estimating  7  on  the  uncertainty  associated  with  the  parameter 
estimate  A.  In  Figure  4.4(a)  we  plotted  the  ratio  SE(A7) / A7  versus  a  for  realizations  of  Dp(cr) 
and  the  ratio  SE(A7)/A7  versus  a  in  Figure  4.4(c).  The  ratios  SE(A7)/A7  in  both  of  these 
examples  appear  to  be  linearly  dependent  on  a.  When  the  full  parameter  set  6  is  estimated, 
we  see  that  the  ratio  SE(A7)/A7  varies  exponentially  with  a  in  Figures  (b)  and  (d)  in  which 
we  plotted  log(SE(d)/d)  versus  a  for  realizations  of  Dp- (a)  and  D™d((j),  respectively. 
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Iq-3  A#  —  Aq  for  Xi  =  0,0.57, 1.14 
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Figure  4.3:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT) ’  0-6  +  jJo >  ■  ■  ■  j  126’  20, 40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57, 1.14.  The  points  denoted 
with  °  are  the  difference  A7  —  Ao,  and  the  points  denoted  *  are  the  difference  A  —  Ao  (a)  The 
result  of  five  realizations  of  D y(cr)  for  values  of  cr  =  0.02,0.05,0.10  (b)  The  result  of  five 
realizations  of  D™nd(a)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 
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Figure  4.4:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  0.6  +  y|g, . . . ,  20, 40, . . . ,  140  and  spatial  nodes  =  0,  0.57, 1.14.  (a)  The  ratio 

SE(A7)/A7  for  five  realizations  of  Djj(cr)  for  values  of  a  =  0.02,0.05,0.10  (b)  The  log  of  the 
ratio  logSE(A)/A  for  five  realizations  of  D.y(cr)  for  values  of  o  =  0.02,  0.05,  0.10.  (c)  The  ra¬ 
tio  SE(A7)/A7  for  five  realizations  of  D)ynd(cr)  for  values  of  a  =  0,  0.02,  0.05,  0.10.  (d)  The  log 
of  the  ratio  log  SE(A)/A  for  five  realizations  of  D-®"nd(cr)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 
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We  also  examined  the  differences  between  7#  —  70  in  Figures  4.5(a)  and  (b).  In  Fig¬ 
ure  4.5(a),  we  plotted  7#  —  70  for  realizations  of  Dy(cr).  In  Figure  4.5(b),  we  plotted  7#  —  70 
for  realizations  of  D™nd(cr).  In  both  Figure  4.5(a)  and  (b),  |7a  —  7o|  is  several  orders  of 
magnitude  less  than  |'y  —  7o|-  Thus  by  estimating  the  entire  parameter  set  9,  we  gain  inac¬ 
curacy  of  our  estimate  of  7.  Moreover,  the  differences  (7  —  7o|  are  an  order  of  magnitude 
greater  than  the  parameter  itself  =  10“3.  In  Figures  4.6(a)-(b),  we  see  that  the  uncer¬ 
tainty  associated  with  the  parameter  estimate  7#  varies  exponentially  with  a.  Figure  4.6(a) 
and  (c),  we  plotted  log(SE(77)/7A)  for  realizations  of  D ij(cr)  and  D™nd(<r),  respectively.  In 
Figures  4.6(b)  and  (d),  we  see  that  log(SE(7)/7)  is  larger  for  realizations  of  DtJ  (a)  (in  (b)) 
than  for  realizations  of  D™ld(cr)  (in  (d))  though  log(SE(7)/7)  (in  Figures  4.6(b)  and  (d)) 
appears  to  be  larger  than  log(SE(7A)/7A)  (in  Figure  4.6(a)  and  (c)). 
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7#  -  7o  for  Xi  =  0, 0.57, 1.14 
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Figure  4.5:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT) !  0.6  +  ifo  j  •  •  • ,  J20  j  20, 40, . . . ,  140  and  spatial  nodes  =  0,  0.57, 1.14.  The  points  denoted 
with  0  are  the  difference  7a  —  70,  and  the  points  denoted  *  are  the  difference  7  —  70  (a)  The 
result  of  five  realizations  of  D ij(cr)  for  values  of  a  =  0.02,0.05,0.10  (b)  The  result  of  five 
realizations  of  D-|nd(a)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 
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Figure  4.6:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT)’  0-6  +  yfo’  •  ■  • ,  yIo’  20, 40, . . . ,  140  and  spatial  nodes  Xj  =  0,  0.57, 1.14.  (a)  The  log  of  the 
ratio  log(SE(7A)/7A)  for  five  realizations  of  D .y(cr)  for  values  of  a  =  0.02,0.05,0.10  (b)  The 
log  of  the  ratio  log  (SE(7)/7)  for  five  realizations  of  D ij(cr)  for  values  of  a  =  0.02,  0.05,  0.10. 
(c)  The  log  of  the  ratio  log(SE(7A)/7A)  for  five  realizations  of  D^nd(cr)  for  values  of  a  = 
0,0.02,0.05,0.10.  (d)  The  log  of  the  ratio  logSE(7)/7  for  five  realizations  of  D-|nd(a)  for 
values  of  a  =  0,  0.02,  0.05,  0.10. 
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We  also  considered  realizations  of  D^cr)  (given  by  (4.9))  and  D'(ind(cr)  (given  by  (4.10)) 
for  spatial  nodes  Xi  =  0,  0.57  and  the  same  temporal  nodes  t3  =  0.6  + 0.6  +  ^E, . . . ,  20, 

40, . . . ,  140.  The  results  of  these  realizations  are  depicted  in  Figures  4.7-4.12. 

We  see  that  for  realizations  of  D,:j(cr)  and  D)jnd(<r)  the  accuracy  of  the  GLS  parameter 
estimation  procedure  for  the  parameter  a  is  similar  for  the  two  random  processes.  We  see 
this  in  Figures  4.7(a)  and  (b)  though  it  does  appear  that  one  realization  of  D-“ld(0.10) 
produced  large  values  of  dA,  &7,  and  a  in  Figure  4.7(b).  In  Figures  4.8(b),  we  see  this 
extreme  realization  of  D+nd(0.10)  also  produced  large  ratios  SE(a#)/a#.  Other  than  this 
extreme  realization,  the  ratios  SE (&#)/&#  for  realizations  of  Djj(cr)  (in  Figure  4.8(a))  and 
D™nd  (cr)  (in  Figure  4.8(b))  appear  to  have  similar  linear  dependence  on  a. 

The  differences  A#  —  A0  are  depicted  in  Figures  4.9(a)  and  (b)  for  realizations  of  D ij(cr) 
and  D)|nd(<r),  respectively.  Again,  we  see  similar  results  for  realizations  of  D; j(a)  and 
realizations  of  D)“ld(<r).  We  also  note,  that  as  was  the  case  for  spatial  nodes  Xi  =  0,  0.57, 1.14 
in  Figure  4.3(a)  and  (b),  |A  —  A0|  is  much  larger  than  |A7  —  A0|  for  both  realizations  of  Dy(er) 
and  realizations  of  D)jnd((j).  The  ratio  SE(A7)/A7  appears  to  vary  linearly  with  a  for  both 
realizations  of  Dy(fj)  (in  Figure  4.10(a))  and  realizations  of  D-“ld(<r)  (in  Figure  4.10(c)), 
while  the  dependence  of  SE(A)/A  is  less  clear  though  it  does  vary  greatly  with  a  for  both 
realizations  D,7  (a)  (log(SE(A)/A)  is  plotted  in  Figure  4.10(b))  and  realizations  of  D((ind(cr) 
(log(SE(A)/A)  is  plotted  in  Figure  4.10(d)). 
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Figure  4.7:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT) !  0.6  +  ifo  ’ ' ' '  ’  lio  ’  40, . . . ,  140  and  spatial  nodes  =  0,  0.57.  The  points  denoted  with 

°  are  the  difference  d7  —  ao,  the  points  denoted  are  the  a\  —  ao,  the  points  denoted  *  are  the 
difference  a  —  ao  (a)  The  result  of  five  realizations  of  Dy(cr)  for  values  of  a  =  0.02,  0.05,  0.10 
(b)  The  result  of  five  realizations  of  D^nd(a)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 
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SE(a#)/a#  for  Xi  =  0,0.57 
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Figure  4.8:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  = 
0.6  +  0.6  +  y§q,  -  -  - ,  1205  20, 40, . . . ,  140  and  spatial  nodes  x,  =  0,  0.57.  The  points  de¬ 

noted  with  0  are  the  ratio  SE(d7)/d7,  the  points  denoted  are  the  SE(dA)/dA,  the  points 
denoted  *  are  the  difference  SE(a)/d  (a)  The  result  of  five  realizations  of  D^(cr)  for  val¬ 
ues  of  a  =  0.02,0.05,0.10  (b)  The  result  of  five  realizations  of  D^nd(cr)  for  values  of 
a  =  0,0.02,0.05,0.10. 
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Figure  4.9:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
IT)  ’  0-6  +  lio  >  ■  ■  ■  >  i2o>  20, 40, . . . ,  140  and  spatial  nodes  =  0,  0.57.  The  points  denoted  with 
°  are  the  difference  A7  —  Ao,  and  the  points  denoted  *  are  the  difference  A  —  Ao  (a)  The 
result  of  five  realizations  of  Djj(cr)  for  values  of  cr  =  0.02,0.05,0.10  (b)  The  result  of  five 
realizations  of  D™nd(a)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 
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Figure  4.10:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
-jTj,  0.6  +  . . . ,  y|q,  20, 40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57.  (a)  The  ratio  SE(A7)/A7 

for  five  realizations  of  Dy(cr)  for  values  of  cr  =  0.02,0.05,0.10  (b)  The  log  of  the  ratio 
logSE(A)/A  for  five  realizations  of  D,:? (a)  for  values  of  a  =  0.02,0.05,0.10.  (c)  The  ratio 
SE(A7)/A7  for  five  realizations  of  D^nd(cr)  for  values  of  a  =  0,0.02,0.05,0.10.  (d)  The  log 
of  the  ratio  log  SE(A)/A  for  five  realizations  of  D-®'nd(cr)  for  values  of  cr  =  0,  0.02,  0.05,  0.10. 
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The  differences  7#  —  70  are  plotted  for  realizations  of  D.,j(cr)  and  D-“ld(<r)  in  Fig¬ 
ures  4.11(a)  and  (b),  respectively.  In  both  figures,  the  values  I7  —  70 1  are  much  larger 
than  the  values  (7^  —  7o|-  This  demonstrates  that  in  these  cases,  estimating  the  full  param¬ 
eter  set  6  contributes  to  inaccuracy  of  the  paramter  estimate  of  7.  Moreover,  the  differences 
I7  —  7o|  are  several  orders  of  magnitude  larger  than  the  “true”  parameter  value  70  =  10~3. 
The  values  of  the  ratios  SE(7#)/7#  vary  a  lot  with  a  in  Figures  4.12(a)-(d).  Figure  4.12(a) 
depicts  log  (SE(7A)/7A)  for  realizations  of  D y(er)  while  D)|nd((T)  depicts  log  (SE(7A)/7A)  for 
realizations  of  D)“ld(a).  In  both  Figure  4.12(a)  and  Figure  4.12(c),  the  ratio  SE(7A)/yA 
appears  to  depend  exponentially  on  cr,  though  the  value  of  SE(7A)/7A  remains  below  one  for 
realizations  that  we  considered.  In  Figures  4.12(b)  and  (d),  we  see  that  log(SE(7)/7)  varies 
between  1-5  for  realizations  of  Dy-  (a)  (in  (b))  and  between  -2-8  for  realizations  of  D™nd(a) 
(in  (d)). 
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Figure  4.11:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
i5o  ’  0-6  +  ifo  >  ■  ■  ■  j  lio  ’  40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57.  The  points  denoted  with 

°  are  the  difference  7a  —  7o,  and  the  points  denoted  *  are  the  difference  7  —  70  (a)  The  result  of 
five  realizations  of  D ij(a)  for  values  of  a  —  0.02,  0.05,  0.10  (b)  The  result  of  five  realizations 


of  D™nd(cr)  for  values  of  a  =  0,  0.02,  0.05,  0.10. 


41 
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Figure  4.12:  The  results  of  realizations  of  simulated  data  with  temporal  nodes  tj  =  0.6  + 
i^o’  0-6  +  y§o  ’  ■  ■  ■  i  i2o>  20,  40, . . . ,  140  and  spatial  nodes  Xi  =  0,  0.57.  (a)  The  ratio  SE(7A)/7>, 
for  five  realizations  of  Dy(cr)  for  values  of  cr  =  0.02,0.05,0.10  (b)  The  log  of  the  ratio 
logSE(7)/7  for  five  realizations  of  Dy(cr)  for  values  of  cr  =  0.02,0.05,0.10.  (c)  The  log  of 
the  ratio  log(SE(7A)/7A)  for  five  realizations  of  D™nd(cr)  for  values  of  cr  =  0,0.02,0.05,0.10. 
(d)  The  log  of  the  ratio  logSE(7)/7  for  five  realizations  of  D™nd(cr)  for  values  of  cr  = 
0,0.02,0.05,0.10. 
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5  Conclusions 


We  investigated  in  [2]  the  behavior  of  solutions  of  (2.2)  (a  partial  differential  equation  on 
a  homogeneous  domain  which  is  derived  from  homogenization  theory)  as  approximations  of 
solutions  of  (2.1)  (a  partial  differential  equation  on  a  randomly  perforated  domain)  in  the 
forward  problem.  Here,  we  considered  the  behavior  of  solutions  of  (2.2)  as  approximations 
of  solutions  of  (2.1)  in  associated  inverse  problems.  Because  it  is  often  difficult  to  verify  a 
priori  the  nature  of  random  error  (whether  or  not  the  error  satisfies  the  OLS  assumptions 
in  (3.3)  or  the  the  GLS  assumptions  in  (4.2)),  we  compared  the  efficacy  of  using  solutions 
of  (2.2)  as  a  model  solution  for  simulated  data  generated  using  solutions  of  (2.2)  to  the 
efficacy  of  using  (2.2)  as  a  model  solution  for  simulated  data  generated  using  solutions 
of  (2.1).  The  results  were  especially  encouraging  for  the  important  parameter  a  (thermal 
diffusivity)  which  will  be  critical  in  our  development  of  NDE  methodology.  The  accuracy  and 
uncertainty  associated  with  the  estimate  of  the  parameter  a  was  similar  for  data  generated 
using  (2.2)  and  (2.1)  for  data  with  relative  and  absolute  added  noise,  for  both  sets  of  spatial 
nodes  Xi  =  0,  0.57  and  Xi  =  0,  0.57, 1.14,  and  for  parameter  sets  6  =  (7,  a,  A),  d7  =  (a7,  A7) 
and  6X  =  (7 Though  estimating  7  presents  difficulties  in  the  inverse  problem  (adds 
uncertainty  and  inaccuracy  to  the  estimate  of  A),  this  affect  was  similar  for  data  generated 
using  solutions  of  (2.2)  and  data  generated  using  solutions  of  (2.1).  There  was  only  one 
example  in  which  there  was  a  significant  difference  between  using  data  generated  using  (2.2) 
and  data  generated  using  (2.1).  For  data  with  added  absolute  random  error  simulated  at 
spatial  nodes  Xi  =  0,0.57,1.14,  the  uncertainty  associated  with  the  parameter  estimate 
was  significantly  larger  for  data  generated  using  (2.1)  than  for  data  generated  using  (2.2). 
We  believe  the  inverse  problem  findings  in  this  report  offer  significant  support  that  such 
methodologies  as  considered  here  will  be  most  useful  in  development  of  NDE  techniques  for 
porous  media  structures. 
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A  Sensitivity 


We  use  the  finite  element  method  to  numerically  solve  (2.1)  and  (2.2).  Here,  we  will  discuss 
the  numerical  solution  of  (2.2)  and  the  sensitivity  functions  in  (3.7);  see  [1]  for  discussion  of 
the  numerical  solution  of  (2.1).  The  finite  element  method  approximates  the  infinite  dimen¬ 
sional  solution  of  a  partial  differential  equation  with  a  finite  dimensional  approximation.  The 
domain  (fl)  is  discretized  using  the  Delaunay  triangulation.  The  finite  dimensional  solution 
is  taken  from  the  space  of  piecewise  two  dimensional  affine  functions,  where  the  solution  is 
affine  on  each  mesh  element  (see  [2],  [14]  and  [17]  for  details).  Specifically,  in  [2],  we  discussed 
the  numerical  approximation  of  U,  the  solution  of  (2.2),  given  by  uN(t,  x)  =  Tj{t)(f)j{x ) 

where  4>j(x)  are  piecewise  affine  basis  element  and  Tj(t)  are  their  time  dependent  coefficients. 
The  coefficients  Tj(t)  are  found  by  solving  the  ordinary  differential  equation  for  T{t)  with 
entries  T)(f) 

PvMft  T(t)  +  (aK  +  XD+  pvlM)  T(t )  =  SfX[toM (t)f,  (A.l) 

where  M  is  an  N  x  N  positive  definite  matrix  with  elements  mt]  =  (0*,  <ff),  K  is  an  N  x  N 
positive  definite  matrix  with  elements  ktl  =  (V0j,  A°V0j),  D  is  an  N  x  N  matrix  with 
components  facfrj  ds ,  /  is  an  Wvector  with  components  /j  =  ,/W4  fiix,  Q)dx  and 

T  is  an  N  column  vector.  To  approximate  Uij{6 *)  in  (3.1),  we  explicitly  integrate  the 
approximation  uN(tj,x )  which  is  piecewise  affine  on  the  source  boundary  uq  so  we  use 

Uij{0*)  «  y  f  l+  uN(tj,s,  0;  9#)  ds. 


Recall  that  in  order  to  calculate  the  covariance  matrices,  we  calculate  the  sensitivity 
matrix 


Xi+m(j—l),k  (d) 


_d_ 

d(  k 


Uij{  C) 


C=e 


Throughout  both  the  generalized  least  squares  and  ordinary  least  squares  parameter  esti¬ 
mation  procedures,  it  is  tacitly  assumed  that  we  use  numerical  estimations  with  reasonable 
convergence  and  that  the  admissible  set  of  parameter  is  compact  and  finite  dimensional  (see 
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[4]  and  [6]).  Given  these  assumptions,  in  both  (3.9)  and  (4.7),  we  may  estimate  x  with  the 
derivative  of  the  numerical  solution  itself.  Explicitly,  we  use 


d  ( 1  [Xi+i 

Xi+m(j-l),k(@)  ~  Xi+m(j—l),k  (^)  =  y~J  J  U  (^ ,  S,  0;  £)  G?S 


(A.2) 


C =o 


Noting  that  the  spatial  nodes  Xi  and  the  interval  width  l  are  parameter  independent,  we 
may  move  the  derivative  under  the  integral  and  replace  uN  with  its  definition 
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Now,  recalling  that  the  basis  elements  <j)j  are  independent  of  6  and  only  Tj(t;  6)  are  dependent 
on  9  =  (7,  a ,  A)  in  (A.l),  we  have 
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d  d  d 

We  need  to  solve  for  —  T7(f;C),  — — T?(f;Q,  and  —Tj(t:C).  First,  in  order  to  calculate 

<77  da  o  A 

d  ->  d  -> 

—Tj(t]9),  we  differentiate  (A.l)  with  respect  to  7.  Let  TJt)  denote  —Tit),  and  recall 

<77  <77 

that  M,  K ,  74,  f  are  independent  of  7.  By  the  chain  rule,  ^( pv^MT(t ))  =  pv'yMT^t)  + 

PvMT{t ),  so  differentiating  (A.l)  we  obtain 


pvM  -^T7(f)  +  (a/(  +  AH  +  py7M)T7(t)  +  pvMT(t)  =  0atxi- 


(A.3) 


In  (A.3)  above,  M,  K  and  D  are  as  defined  in  (A.l),  and  0atxi  is  the  iV  x  1  vector  of  zeros. 
The  ordinary  differential  equation  (A.3)  has  a  term  that  involves  T(f),  so  (A.3)  and  (A.l) 
must  be  solved  simultaneously.  Similarly,  we  take  the  derivative  of  (A.l)  with  respect  to  a 
to  obtain 

pvM  ^ Ta(t )  +  (, aK  +  A  D+  pvjM)fa(t )  +  KT{t)  =  O^xi,  (A.4) 

which  must  also  be  solved  with  (A.l).  Finally,  we  take  the  derivative  of  (A.l)  with  respect 
to  A  to  find  the  system  of  ordinary  differential  equations  for  T\  given  by 

pvM  ^Tx(t)  +  (aK  +  A  D  +  pvjM)fx(t)  +  DT(t)  =  0jvxl,  (A. 5) 

which  also  must  be  solved  simultaneously  with  (A.l). 
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